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Abstract. The stochastic mutual repressor model is analysed using perturba- 
tion methods. This simple model of a gene circuit consists of two genes and three 
promoter states. Either of the two protein products can dimerize, forming a re- 
pressor molecule that binds to the promotor of the other gene. When the repressor 
is bound to a promotor, the corresponding gene is not transcribed and no protein 
is produced. Either one of the promoters can be repressed at any given time or 
both can be unrepressed, leaving three possible promotor states. This model is 
analysed in its bistable regime in which the deterministic limit exhibits two stable 
fixed points and an unstable saddle, and the case of small noise is considered. On 
small time scales, the stochastic process fluctuates near one of the stable fixed 
points, and on large time scales, a metastable transition can occur, where fluctua- 
tions drive the system past the unstable saddle to the other stable fixed point. To 
explore how different intrinsic noise sources affect these transitions, fluctuations 
in protein production and degradation are eliminated, leaving fluctuations in the 
promotor state as the only source of noise in the system. Perturbation methods 
are then used to compute the stability landscape and the distribution of transi- 
tion times, or first exit time density. To understand how protein noise affects the 
system, small magnitude fluctuations are added back into the process, and the 
stability landscape is compared to that of the process without protein noise. It is 
found that signiflcant differences in the random process emerge in the presence of 
protein noise. 



1. Introduction 

Random molecular interactions can have profound effects on gene expression. Because 
the expression of a gene can be regulated by a single promotor, and because the 
number of mRNA copies and protein molecules is often small, deterministic models 
of gene expression can miss important behaviors. A deterministic model might show 
multiple possible stable behaviors, any of which can be realized depending on the initial 
conditions of the system. Different stable behavior that depend on initial conditions 
allows for variability in response and adaptation to environmental conditions [1] . 

Although in some cases, noise from multiple sources can push the behavior far 
from the deterministic model, here we focus on situation where the system fluctuates 
close to the deterministic trajectory (i.e., weak noise). Of particular interest is 
behavior predicted by a stochastic model that is qualitatively different from its 
deterministic counterpart [2] , even if the fluctuations are small. 

Several interesting questions emerge when including stochastic effects in a model 
of gene expression. For example, what are the different sources of fluctuations affecting 
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a gene circuit? Can noise be harnessed for useful purpose, and if so, what new functions 
can noise bring to the gene-regulation toolbox? One way in which noise can induce 
qualitatively different behavior occurs when a rare sequence of random events pushes 
the system far enough away from one of the stable deterministic behaviors that the 
system transitions toward a different stable dynamic behavior, one that would never 
be realized in the deterministic model without changing the initial conditions. For 
example, if the deterministic model is bistable, fluctuations can cause the protein 
concentration to shift between the different metastable protein concentrations. This 
happens when fluctuations push the system past the unstable fixed point that separates 
two stable fixed points. 

While often times a spontaneous change in gene expression might be harmfull, it 
might also be beneficial. For example, in certain types of bacteria, a few individuals 
within a population enter a slow-growth state in order to resist exposure to antibiotics. 
In a developing organism, a population of differentiating cells might first randomly 
choose between two or more expression profiles during their development and then later 
segregate into distinct groups by chemotaxis. In both examples, switching between 
metastable states leads to mixed populations of phenotypic expression [3]. This leads 
to the question of how cells coordinate and regulate different sources of biochemical 
fluctuations, or noise, to function within a genetic circuit. 

In many cases, the genes within a given circuit are turned on and off by regulator 
proteins, which are often the gene products of the circuit. If a gene is switched on, its 
DNA is transcribed into one or more mRNA copies, which are in turn translated 
into large numbers of proteins. Typically, the protein products form complexes 
with each other or with other proteins that bind to regulatory DNA sequences, or 
operators, to alter the expression state of a gene. For example, a repressor binds to an 
operator which blocks the promotor — the region of DNA that a polymerase protein 
binds to before transcribing the gene — so that the gene is turned off and no mRNA 
are transcribed. This feedback enables a cell to regulate gene expression, and often 
multiple genes interact within groups to form gene circuits. 

Understanding how different noise sources affect the behavior of a gene circuit and 
comparing this with how the circuit behaves with multiple noise sources is essential for 
understanding how a cell can use different sources of noise productively. Fluctuations 
arising from the biochemical reactions involving the DNA, mRNA, and proteins are 
commonly classified as "intrinsic" noise [4] . One important source of intrinsic noise is 
fluctuations from mRNA transcription, protein translation, and degradation of both 
mRNA and protein product. This type of noise is common among many of the 
biochemical reactions within a cell, and its effect is reduced as the number of reacting 
species within a given volume grows large. Another source of intrinsic noise is in the 
expression state of the genes within the circuit. Typically there is only one or two 
copies of a gene within a cell, which means that thermal fluctuations within reactions 
with regulatory proteins have a significant effect on mRNA production. Here, we 
consider the situation where transitions in the behavior of a gene circuit are primarily 
driven by fluctuations in the on/off state of its promotor and examine the effect of 
removing all other sources of noise. 

Stochastic gene circuits are typically modelled using a discrete Markov process, 
which tracks the random number of mRNA and/or proteins along with the state of 
one or more promotors [5-9] (but see also [10]). Monte-Carlo simulations using the 
Gillespie algorithm can be used to generate exact realizations of the random process. 
The process can also be described by its probability density function, which satisfies 
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a system of linear ordinary differential equations known as the Master equation. The 
dimension of the Master equation is the number of possible states the system can 
occupy, which can be quite large, leading to the problem of dimensionality when 
analyzing the Master equation directly. However, for the problem considered here, the 
full solution to the Master equation is not necessary in order to understand metastable 
transitions. 

The motivating biological question we consider here is what percentage of a 
population of cells can be expected to exhibit a metastable transition within a given 
timeframe. If a spontaneous transition is harmfuU to the cell, one expects that reaction 
rates and protein/DNA interactions should evolve so that transition times are likely 
to be much larger than the lifetime of the cell. On the other hand, if spontaneous 
transition are functional, transition times should be tuned to achieve the desired 
population in which the transition occurs. In either case, the key quantity of interest is 
the distribution of transition times between metastable states, regardless of the noise 
source driving the transition. 

Except for a few special cases [9, 11] exact results, even for the mean transition 
time, are not possible and approximation techniques or Monte-Carlo simulations must 
be used. However, because rare events typically involve long simulation times where 
large numbers of jumps occur, Monte-Carlo simulations are computationally expensive 
to perform, leaving perturbation analysis ideally suited for the task. 

Past studies of metastable transitions, where perturbation methods are applied 
to the Master equation, have used a simplifying assumption so that the state of the 
promotor is not accounted for explicitly [6,7]. The assumption is that proteins are 
produced in "bursts" during which one or more mRNA copies are translated to rapidly 
produce many proteins. In these models, production bursts occur as instantaneous 
jumps, with a predefined distribution determining the number of proteins produced 
during a given burst. More recently, Assaf and coworkers analyzed a model where 
the on/off state of a single stochastic promotor is accounted for explicitly, and 
mRNA copies are produced stochastically at a certain rate when the promotor is 
turned on [12]. However, the case where the model contains an arbitrary number 
of promotors or promotor states has not been addressed, and as we show in this 
paper, accounting for even just three promotor states is nontrivial. Similar asymptotic 
methods have also been developed to study metastable transitions in continuous 
Markov processes [13-18], but cannot be applied to a discrete chemical reaction 
system because continuous approximations, such as the system-size expansion, of 
discrete Markov processes do not, in general, accurately capture transition times [19]. 
Another source of difficulty that arises from isolating promotor state fluctuations as 
the only source of noise is that the resulting state space of the Markov process is both 
continuous and discrete. 

After removing all sources of intrinsic noise except for the fluctuating promotor — 
by taking the thermodynamic limit — the protein levels change deterministically and 
continuously, and the promotor's state jumps at exponentially-distributed random 
times. The random jumps in the promotor's state makes the protein levels appear 
random, even though they are only responding deterministically to changes in the 
promotor state. Such random processes are sometimes called hybrid systems, or 
piecewise deterministic [20,21]. Here we refer to it as the quasi-deterministic (QD) 
process because we are taking part of the randomly fluctuating discrete state of the 
system (the number of protein molecules) and replacing it with a deterministically- 
changing continuous state. Recently, we have developed asymptotic methods for 
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metastable transitions — similar to those applied to the discrete Master equation — for 
Markov processes with both discrete and continuous state spaces [22,23]. However, 
these methods do not account for two or more continuous state variables, which 
restricts the genetic circuit the method can analyze to one with a single protein 
product. 

In this paper, we develop new perturbation methods so that we can study 
metastable transitions in genetic circuits driven by promotor fluctuations. These 
methods are based on previous theory developed for one-dimensional velocity 
jump processes, and are generalized to account for the multiple continuous states 
representing the quantity of proteins produced by the genetic circuit. They also fit 
within a larger framework of methods to study metastable transitions in continuous 
Markov processes [18] and in discrete Markov processes [6,7,12,24-30]. For illustration, 
we use a simple model known as the "mutual repressor" model [31], which contains 
two genes, two promotors, and three promotor states. Although our example 
considers only three promotor states, the methods presented are general and can 
account for an arbitrary number of promotor states. For a range of parameter 
values, the deterministic limit of the mutual repressor model is bistable, having 
two stable fixed points separated by an unstable saddle point. For the stochastic 
model, the deterministic forces create two confining wells surrounding each stable 
fixed point, separated by a stability barrier along the separatrix that contains the 
unstable saddle. This geometric interpretation is given by taking the logarithm of the 
stationary probability density function, which we refer as the "stability landscape." 
An approximation of the first exit time density is found for the random process to 
escape over the stability barrier from one of the stability wells to the other. Using 
this model, we seek to answer the following question. How does the random process 
change when protein noise is removed, leaving the state of the promotor as the only 
source of randomness? That is, are there any qualitative differences in the behavior 
of the system without other sources of intrinsic noise? 

The paper is organized as follows. In Section 2 the Mutual Repressor model 
is presented along with its reduction to the QD process, and then in Section 3, 
perturbation methods for estimating the exit time density are applied to the QD 
process. For comparison, the stability landscape is also computed for the full process in 
Section 3.1.2, which includes fluctuations in protein production/degradation. Finally, 
results are presented in Section 4, and the QD process is compared to the full process, 
using analytical/numerical approximations and Monte-Carlo simulations. 

2. Mutual repressor model 

The mutual repressor model [31] is a hypothetical gene circuit consisting of a single 
promotor driving the expression of two genes: N and M. Each protein product can 
dimerize and bind to the promotor to repress the expression of the other. When no 
dimer is bound to the promotor, both genes are expressed equally. Thus, the promotor 
can be in one of three states: bound to a dimer of protein one Oq, unbound Oi, or 
bound to a dimer of protein two O2 ■ Let the number of protein product of gene N and 
M be n and m, respectively. It is assumed that the mRNA and protein production 
steps can be combined into a single protein production rate and that the dimerization 
reaction is fast so that it can be taken to be in quasi-steady-state, we then have the 
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following transition between the three promotor states 

Oo ^ Oi O2, (2.1) 

where k is a rate and /3 is a nondimensional dissociation constant. Protein TV (M) is 
produced at a rate a while the promotor is in states Oq,! ((^1,2), and both proteins 
are degraded at a rate 6 in all three promotor states. 

The probability density function pk{n,m,t), for promotor state fc = 0,1,2 and 
protein numbers n, m = 0, 1, 2, • • •, satisfies the Master equation 

^p{n,m,t) = [A + W]p, (2.2) 

where 

" -/3 n{n~ 1) 

A = K /3 -n{n ~ 1) - m{m - 1) fi (2.3) 

m(m — 1) —(3 

is the matrix responsible for promotor state transitions. The diagonal matrix W, 
responsible for changes in protein numbers, has elements 

Wo =D + a(E; - 1), (2.4) 

Wi==D + a(E-+E„-2), (2.5) 

W2-D + a(E;,-l), (2.6) 

with 

D = 5 [(E+ - l)n + (E+ - l)m] . (2.7) 
The shift operators E^ are defined according to 

%*/(j) =/(j±l)- (2.8) 

We now introduce the nondimensional variables t ^ tS, n ^ x/j, and m — y/j, 
where 'j = S/a. Then, the Master equation (2.2) for the rescaled probability density, 
p{n,m,t) — )■ p{x,y,t), becomes 



dt 



p{x,y,t) 



1 



1. 



-W 



(2.9) 



with dimensionless parameters b 

-b 



and 6 = ^ 



Ay — 



x{x + 7) 
b —x{x + j) — y{y + -f) 

y{y + i) 



The matrices are given by 


b 



and 



M/ = ((e^^ - l)a; + (e^^ - l)y)/ 



dy 



(2.10) 



(2.11) 



diag(e' 



~dx 



l,e 



-dx 



+ e~^y -2,e-'^y - I). 



The operators e^^^ and e*^^ are defined in terms of Taylor series expansions (in small 
7) which replace the shift operators E^ „j with 



(±7)" 9"/ 
n! dx^ 



/(x±7). 



(2.12) 
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Assume that 7 ^ 1 is a small parameter, so that there is a large average 
number of proteins, and assume also that the parameter e <C 1 is small, which 
reflects rapid switching between promotor states compared to the rate of protein 
production/degradation. Because we have two small parameters in our system, 7 
and e, when perusing an asymptotic solution, we must carefully consider how the 
limit e — > and 7 — > is taken, or more practically, how large e is compared to 7. 
The fluctuations in the promotor state are controlled by e, and in the limit e — ?► 0, 
the transitions are infinitely fast so that the promotor behaves deterministically. The 
fluctuations in protein levels is controlled by 7, and in the limit 7 — >■ the protein 
production/degradation behaves deterministically. Since we are concerned primarily 
with rare transitions driven by promotor fluctuations and not by fluctuations in the 
protein production/degradation reaction, we assume that 7^6 (i.e., ^ ^ 1). 

Taking both limits, e — )■ and 7 — 0, yields the fully-deterministic dynamics, 

i = fix, y), y = f{y, x), (2.13) 

where 

f{x,y)^ ^ \. -X. (2.14) 

Note the symmetry in the problem; the deterministic system is unchanged if we 
exchange x ^ y. Dynamically, the system is bistable for < 6 < At 5 = 6„ = 4/9 
there is a saddle-node bifurcation, and for b > bu there is a single stable fixed point, 
we consider the case of bistability and chose b <^ bu- In Fig. 1 the nullclines and fixed 
points are shown. The two stable fixed points are located near the corners, and the 



y 




^'8.0 0.2 0.4 0.6 0.8 1.0 

X 

Figure 1. Deterministic dynamics for b = 0.15. The black curve shows the 
y-nuUcline and the grey curve shows the x-nullchne. The green circles show the 
stable fixed points, the red circle shows the unstable saddle. The blue curve 
shows a stochastic trajectory leaving the lower basin of attraction to reach the 
separatrix. 

unstable saddle point is located along the separatrix. Arrows show the eigenvectors 
of the Jacobian with their direction determined by the sign of the eigenvalues, all 
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of which are real. A stochastic trajectory that starts at the lower stable fixed point 
remains nearby for a long period of time until a rare sequence of jumps carries it to 
the separatrix. Because the separatrix is the stable manifold, trajectories are most 
likely to exit near the unstable saddle point. 

To remove protein noise from the system, consider the limit 7—^0, with e ^ 1 
fixed, so that the protein production/degradation process is deterministic within each 
promotor state, while the promotor state remains random. The master equation (2.9) 
becomes 



dp 

dt 



d_ 
dx 



{Fp) 



d_ 
dy 



{Gp) + -Ap, (x,y)e(0,l) 



(2.15) 



where 



and 



F{x, y) = diag(l-a;, 1-x, -x), G{x, y) = diag(-?/, 1-y, l-?/).(2.16) 



(2.17) 





' -b 


x2 







h 




b 









-b 



3. Transition rate theory 

The focus of the remaining analysis is to obtain an accurate approximation of the first 
exit time density function (FETD) for the QD process to evolve from one metastable 
state to the another. To obtain the FETD, we supplement an absorbing boundary 
condition to the governing equation along the separatrix, 

r={ix,y):0<x = y<l}, (3.1) 

of the deterministic dynamics, which is the barrier the process must surmount in 
order to transition to the other metastable state. For the Master equation (2.2), the 
absorbing boundary condition is simply 

p(n,n,t) = 0, for n==0, 1,---. (3.2) 

For the QD CK equation (2.15), the absorbing boundary condition is 

Pn{x,y)^0, for (x,?/) e r, and (i^„^„,G'„,„) • n < 0, (3.3) 

where n is the unit vector normal to the boundary F. Then, the domain for the QD 
process with the absorbing boundary is given by 

X> EE (0, 1)2 n {x < y}. (3.4) 

Note the choice of the lower triangular region (instead of the upper triangular region 
with X > y) is arbitrary due to the symmetry in the problem. 

To see how the absorbing boundary on F sets up the exit time problem, define 
T to be the random time at which the separatrix is reached for the first time, given 
that the process starts at the stable fixed point (x*,?/*) e V. Consider the survival 
probability 

3 

S{t) = Yl / Pn{x,y,t)dA, (3.5) 
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which is the probability that t < T. The FETD (or probability density function for 
T) is then 

m - -% (3^6) 

The FETD for the QD process can be approximated using perturbation methods 
as follows. Suppose we have a CK equation of the form 

dtp = -C,p, (3.7) 

where is a linear operator acting on the continuous and discrete state variables of 
the density function. In the case of the QD CK equation (2.15), we have 

C,^—{Fp) + —{Gp)^-Ap, ix,y)&V (3.8) 

Assume that £^ has a complete set of eigenfunctions, {<pj}, so that the solution can 
be written 



oo 



p(x,y,i) =^Cj0,(a;,y)e-^^*, (3.9) 

for some constants Cj, and that all of the eigenvalues, Aj, are nonnegative. Assume 
further that if we impose a reflecting boundary condition on F then the principal 
eigenvalue, Aq = 0, is the only zero eigenvalue and the eigenfunction, (pQ, is the 
stationary density of the process (after appropriate normalization). Furthermore, we 
assume that the stationary density is exponentially small on the boundary. Then, 
if instead we place an absorbing boundary condition on F, the stationary density 
no longer exists and the principal eigenvalue is exponentially small in e, while the 
remaining eigenvalues are much larger so that the solution resembles the stationary 
density after some initial transients. It is this difference in time scales that we exploit 
to approximate the FETD. 

A universal feature of the FETD for rare events is its exponential form (which 
follows from the separation of time scales) because the time dependence is e^'*'"*, for 
Xit > 1. Indeed the FETD (3.6) is 

T{t) Xae'^"\ forAiOl. (3.10) 

Thus, for large times the FETD is completely characterized by the principal eigenvalue, 
Aq. The mean exit time is simply 1/Ao, which means that the eigenvalue also has the 
physical interpretation of the rate at which metastable transitions occur. To obtain an 
approximation of this eigenvalue, we use a spectral projection method, which makes 
use of the adjoint operator C*. Consider the adjoint eigenfunctions, {^j}, j = 0, 1, • • •, 
satisfying 

C:^, = AC,, (3.11) 

and (0i,Cj) — so that the two sets of eigenfunctions are biorthogonal. 

Soppose that the boundary condition is ignored and <po is approximated by the 
stationary density, p. By application of the divergence theorem, the adjoint operator 
is such that 

(p, C:^o) = {CeP, 6) + £ ^oiF - G)p ds (3.12) 
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where the boundary contribution is nonzero because p does not satisfy the absorbing 
boundary condition. Then, since = 0, the principal eigenvalue is 

In the remainder of this section, we use (3.13) to approximate the principle eigenvalue 
by approximating the stationary density, p, and the adjoint eigenfunction, ^q- 

3.1. Quasi- stationary distribution 

In this section, we obtain an approximation to the quasi-stationary density, p{x,y), 
using a WKB approximation method, we begin by illustrating the procedure for the 
QD process. That is, we seek an approximation of the solution to the equation 



e ox ay 



Pix,y)=0. (3.14) 



Consider the anzatz 

"1 



Pix,y) = (ro(2^,y) + eri(a;,y))exp 



e 



^ix,y) + k{x,y) 



(3.15) 



where ro,i(a;,?/) are 3-vectors and both <^{x,y) and k{x,y) are scalar functions. Note 
that in other studies of gene regulation models where similar methods are used, the 
small parameter in the exponential is 7. This difference in scaling arrises from the 
assumption that the metastable transitions are driven by fluctuation in the promotor 
and not the production of protein. Substituting (3.15) into (3.14) and collecting 
leading order terms in e yields 

[A+pF + qG]ro = 0, (3.16) 

where 

The prefactor term Tq (up to a normalization factor) is simply the nullspace of the 
matrix M = [A+pF+qG], and we assume that it is normalized so that X]n=o(''o)ji = 1- 
Using Theorem 3.1 in Ref. [23], we can provide necessary and sufficient conditions for 
Tq to be unique and positive. For any fixed {x,y,p, q), there exists a unique vector 
To > 0, satisfying (3.16) if and only if the diagonal matrix H = pF + qG is such that 
at least two of its elements have oposite sign. That is, there exist i, j, with i 7^ j, such 
that Hi^iHjj < 0. It is interesting to speculate that once the solution {p{x, y), q{x, y)) 
to (3.60) is substituted into the matrix M that this requirement is satisfies for all 
{x,y) G (0,00)^. However, this is not necessarily the case, which means that the 
quasi-stationary density is restricted to a subdomain where ro{x,y) > 0. 

It is obvious that the protein levels must be bounded within the domain {x,y) S 
(0,1)^ when protein production/degradation is deterministic. That is, if the gene 
remains in the unrepressed state, the protein level tends toward the mean value 
{n,m = 1/7 or a:, 2/ — 1) but never exceeds it since protein levels do not fluctuate 
(unless the promotor state fluctuates). However, it is not as obvious that the total 
amount of protein is further bounded so that 1 < x + y < 2, which means the domain 
is further restricted to the upper triangular portion of the unit square. This means 
that once a trajectory enters, it remains in this domain for all time and cannot escape. 
To show this, we need simply look at the rate of protein production/degradation for 
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each state normal to the Hne y — 1 — x. This gives the rate for each of the promotor 
states {x on, both on, y on) when both protein levels satisfy y = 1 — x. These rates 
are given by the diagonal components of the matrix F{x) + G(l — x) = diag(0, 1, 0), 
where F and G are defined in (2.16). It is evident that when no repressor is bound 
and both proteins are produced, the flux across this line is in the positive direction, 
and when one repressor is bound, there is no flux across this line. 

The leading order result (3.16) defines a nonlinear partial differential equation for 

^, 

n{x, y,p, q) = det[A{x, y) + pF{x, y) + qG{x, y)] - 0. (3.18) 

The function Ti. is referred to as the Hamiltonian for the system, due to the similarity 
to classical hamiltonian dynamics. An implicit assumption, 

dy dx ^ ^ 

is present to ensure that {p, q) is the gradient of the scalar field $. The above system 
can be solved by the method of characteristics to obtain the system of ordinary 
differential equations (see [32] pg. 360), 

on . dn 

where each variable is parameterized by t, which should not be confused with physical 
time. The above system of ordinary differential equations is supplemented with an 
equation for the stability lanscape 

^ d<i> . . 

dn on 

and solutions specify $ along a curve in the plane {x{t),y(t)). A family of curves is 
defined by specifying Cauchy data 

x{0)=xoi9), yiO)=yo{9), p(0)=Po(^?), g(0) - <Zo(e), (3.24) 

along a curve parameterized by 0. 

One of the difficulties found in this method is determining Cauchy data. At the 
stable fixed points, the value of each of the variables is known (i.e., p ~ q — and 
X = x^,,y — y^,) but data at a single point cannot hope to generate a family of rays. 
Therefore, data must be specified on an ellipse surrounding the fixed point, using the 
Hessian matrix, 

dydx dy^ 

Expanding the function $ in a Taylor series around the fixed point yields the quadratic 
form, 

$(a;,y)« Vzr, r= ^^M, (3.26) 
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as its leading order term. Cauchy data is specified on the ellipse 



lri0)Zri0f 



(3.27) 



for some suitably small uj <^ I. In practice, a; must be small enough to generate 
accurate numerical results, but large enough so that trajectories can be generated to 
cover the domain. On the elliptical contour, the initial values for p, q are 

It can be shown [18] that the Hessian matrix is the solution to the algebraic 
Riccati equation, 

ZBZ + ZC + C'^Z = 0, 



(3.28) 



where 



B = 



dqdp 



dpdq 



dpdx 
dqdx 



dpdy 
dqdy _ 



(3.29) 



(3.30) 



evaluated at p = q = 0, x = x*, and y = y*. This equation can be transformed into a 
linear problem (in order to actually solve it) by making the substitution Q = Z~^ to 
get 

B + CQ + QC^ = 0. (3.31) 



3.1.1. An equation for k{x,y) An equation for the scalar function k{x,y) is found 
by substituting (3.15) into (3.14) and keeping second order terms in e, to get 

[A+pF + qG] ri - £(Fro) + ^(Gro) - i^^F - ^Gj ro.(3.32) 

For solutions ri to exist, the Fredholm Alternative Theorem requires that for 

l^[A + pF + qG]=0 (3.33) 

it must be that 



I' 



'-iFro) + liGro)-('J^F-'J^G 



(3.35) 



dx"^ ' dy^"" \dx dy"" ' ^ ^ ("^ ) 

Note that if rg spans the right nuUspace oi A + pF + gG, then the left nullspace is 
also one dimensional. After rewriting (3.34), we have the PDF for k given by 

Although the solution to this equation can be formulated by the method of 
characteristics, it requires values of the vectors Tq and i, which in turn require the 
solution to the ray equations (3.20). Since rays must be integrated numerically in most 
cases, solving (3.35) along its own characteristics is impractical. However, (3.35) can 
be computed along the characteristic curves of (3.20) as follows. First, differentiating 
(3.35) along characteristics yields 



; dk . 

ox 



dk 
dy 



y 



dkdn 
dx dp 



dkdn 
dy dq ' 



(3.36) 
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Using the fact that {l'^Gro)x — {l'^Fro)y = along characteristics, we can define 



h{x,y) 



du 

dp 



dU 
dq 



Then, after combining (3.35) and (3.36) we have that 



h{x,y)f 



dF 
dx 



dG 
dy 



Vq + F 



dro dro 



dx 



G 



dy 



(3.37) 



(3.38) 



The above requires values of ^ and which are not provided by the system (3.20). 
To obtain these, a formula is needed to relate the Hessian matrix, Z{x,y), of ^{x,y) 
to Vfq. Then, the Hessian matrix can be computed by expanding the system of ray 
equations (3.20). 

First, differentiate both sides of equation (3.16) to get 

[A+pF + qG] Vro = - (VA + W{pF) + V(gG)) ro, (3.39) 

The Fredholm Alternative Theorem requires that 

F {VA + V{pF) + V{qG)) ro - F(VAf)ro = 0, (3.40) 

which is always true since Mtq — M'^l — and 

= V(FMro) (3.41) 

= (V«^)Mro + F(VM)ro + FM{Vro) (3.42) 

= F(VM)ro. (3.43) 

The general solution to (3.39) is 

Vro = -A/1'VMro + OTq, (3.44) 

where is the pseudoinverse of the matrix M and a is an unknown constant. Since 
the vector ro is normalized so that its entries sum to one, it follows that J2ni'^^o)n — 0. 
Summing over both sides of equation (3.44) then yields 



2 

E 

n=0 



(AftvA/ro)„. 



Thus, we have that 



Vrn 



2 

E 



(3.45) 



(3.46) 



Equation (3.46) gives a relationship between Z, Tq, and Vro. To obtain the 
Hessian matrix, Z(x,y), away from the fixed point, the ray equations are extended to 
include the variables 

dy 



Pj 



dx 
duj' 
dp 
dui ' 



1j 



dv 



(3.47) 
(3.48) 



for j = 1,2. A good choice for the current problem is to take ui = xq{9) and 
U2 — ya{0), where {xo{9) , yo{0)) is a point on the initial curve defined by (3.27). The 
Hessian matrix is then obtained using 

(3.49) 



Pi 


P2 


= z 


Xi 


X2 


. 91 


q2 






y2 
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As long as the matrix on the RHS is invertible, the matrix Z can be obtained along 
characteristics and k can be integrated numerically using (3.38) and (3.39). The 
dynamics for the extended variables (3.47) is given by 

^ J{x,y,p,q)vj, (3.50) 

where Vj = {xj,yj,pj,qj)^ and J{x,y) is the Jaciobian matrix for the system (3.20). 
One can choose different variables Uj with which to extend the system based on what 
works in practice, and the only thing one must change arc the initial conditions. For 
our choice, the initial conditions are 

x,{0) = 6,^j, y,(0)=<52,„ (3.51) 

Pj(0) = Zj4(x*,?/*), qj{0) = Zj^2{x*,y*), (3.52) 

where the matrix Z{x^,,y^,) is the solution to (3.29). 



3.1.2. Stability landscape for the full process The above analysis can be repeated to 
obtain a Hamiltonian system for the full process (with protein noise), but a choice 
must be made for how the limit e — >■ 0, 7 — ;> is taken. First consider the equation 
for the quasi-stationary density, p, for the full process (2.9) 



1 , 1 

-A + -T 
e 7 _ 



p{x,y) = 0. 



(3.53) 



Here, the domain is the cone < a; < y < 00. The quasi-stationary density is assumed 
to have the form 

-Hx,y) , (3.54) 



pix,y) = r(x,y)exp 



where again r is a 3-vector and $ is a scalar function representing the stability 
landscape. Note that we have ignored higher order terms here because we only want 
the Hamiltonian function for comparison to the QD process. Substituting (3.54) into 
(3.53) does not lead to any meaningful equation for $ at leading order unless we make 
an assumption about how the limit 7 — ^ is taken. There are two relevant cases: 
7 — o(e) and 7 — 0(e). In the former case, one recovers the QD process (3.16), and 
in the latter case, collecting terms of leading order in e, with 7 = ipe, yields 

^-H r = 0, 



A- 



where H = di&g{Ho{x,y,p,q),Hi{x,y,p,q),H2{x,y,p,q)) and 



Ho{x,y,p,q) 
Hi{x,y,p,q) 
H2{x,y,p,q) 

h{x,y,p, q) 



hix,y,p,q) 

h'{x,y,p,q) 

h{x,y,p,q) 
'pVP - 1 



1 



-vg 



1 



Thus, the Hamiltonian for the full process is 

'H{x,y,p,q) = det [A{x,y) + Ii{x,y,p,q)] 
which we refer to as the full Hamiltonian 



(3.55) 

(3.56) 
(3.57) 
(3.58) 
(3.59) 

(3.60) 
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The differences between the full process and the QD process is nicely illustrated 
by comparing their associated Hamiltonians. Notice that the full Hamiltonian (3.60) 
is a transcendental function of p and q, whereas the Hamiltonian for the QD process 
(3.18) is a cubic polynomial in p and q. One can view this as a Taylor series expansion 
of the full Hamiltonian about (p{p,q) — (0,0). For this reason, the QD process, as 
an approximation for the full process with a small amount of protein noise, is only 
valid within a neighborhood of a deterministic fixed point. This is, of course, just a 
reflection of the fluctuation dissipation theorem [33] . 

An example of numerical integration (for details regarding numerics see the 
Appendix) of the ray equations (3.20) for the QD (3.18) and full (3.60) Hamiltonian is 
shown in Fig. 2. The QD rays are shown above the separatrix for comparison. Notice 



Figure 2. Numerical solutions to the ray equations for the QD and full 
Hamiltonians. The grey lines are characteristic projections, {x{t),y{t)). The 
triangular restricted domain boundary for the QD process is shown with black 
lines. The two symmetric stable fixed points are represented by green circles, and 
the unstable fixed point by a red circle. The separatrix along which an absorbing 
boundary condition is imposed is shown by a dashed black line. Parameter values 
used are 6 = 1.5, oj = 10~^ for the QD Hamiltonian, and u) = 10"^'' for the full 
Hamiltonian. 



that the QD rays are contained within a triangular domain, while the rays from the 
full Hamiltonian cover the entire domain. This is due to the domain restriction that 
occurs when removing the protein fluctuations from the process. 

3.2. Adjoint eigenfunction 

Up to terms that are exponentially small in e, the adjoint eigenfunction, 6, satisfies 



Characteristics 




X 



To make things easier, we change coordinates with 




(3.61) 




(3.62) 



Isolating intrinsic noise sources in a stochastic genetic switch 



15 



so that 

T = x + y — l, a~x~y. (3.63) 

This transforms the absorbing boundary, a; — y = 0, to the vertical Une, a — Q. Then, 
(3.61) becomes 



lA^{r,a) + F{T)-l-^G{a)^ 



^o(T,a) = 0. (3.64) 



where 



A{T,a) = A{x{T,a),y{T,a)), (3.65) 
F(t) = F + G = Amg{-T,l-T,-T), (3.66) 
G(cr) =i^-G = diag(l-(T,-(T,-(l + cr)), (3.67) 

where A, F , and G are defined in (2.16) and (2.17). The absorbing boundary condition 
is then 

^(°)(r,0) = 0, for re (-1,1). (3.68) 

Before proceeding, it is convenient to make the foUowing definitions. In the rest 
of this section, we make frequent use of the eigenvectors (right eigenvectors i/'„ and 
left eigenvectors J7„) and eigenvalues, /i„, satisfying 

i(T,CT)l/?„(T, cr) = ^„(t, cr)G'((T)'0„(T,cr) (3.69) 
A(t, afrin{T, a) = ^„(t, cr)G'(cr)?7„(r, a). (3.70) 

We normalize the two sets of eigenvectors (which are biorthogonal) so that ijfGxpj = 
6ij. Note that because the matrices A and G are functions of (t, ct), so are the 
eigenpairs. It is easily shown that one of the eigenvalues is zero for all values of (r, a), 
which we set to /io = 0. The right eigenvector, if)o, is then given by the nuUspace of 
the matrix A 




Furthermore, the corresponding left eigenvector is simply 

r7o = l= (1,1,1)^. (3.72) 
It is convenient to define distinct notation for the eigenpairs evaluated on F, with 

^„(r) = ■j/)„(t, 0), ■i7„(r) = ?7„(t,0), /i„(T) ee ^„(t, 0). (3.73) 

At the boundary one of the eigenvalues, fii say, vanishes and the eigenspace for the 
zero eigenvalue is degenerate (i.e. there are two zero eigenvalues but the nuUspace is 
one dimensional) which means that fli = 0, xpi — "00 and fji — 1. 

The approximation of the adjoint eigenfunction proceeds using singular 
perturbation methods, along the lines of [22]. Three solutions are found which are 
valid in different regions of the domain: an outer solution, a boundary layer solution 
for the 0(e) strip near the absorbing boundary, and a transition layer solution in the 
0(e^/2) overlap region between the other two. 

Away from the boundary, the exact solution (that does not satisfy the boundary 
condition) is 

lout = 1. (3.74) 
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To obtain a uniform asymptotic approximation that also satisfies (3.68), a boundary- 
layer solution is needed. Consider the stretched variable z = cr/e. To leading order in 
e, the boundary-layer solutions, ^bi(''', z), satisfies 

G(0)^+i(r, 0)^^1 = 0, (3.75) 
oz 

where (5(0) = diag(l, 0, —1). The solution has the form 

|bi(r, z) = col + ci(C + zl) + C2i]2e-~^'', (3.76) 

where 172, fj.2 is the only eigenpair (on the boundary) with a nonzero eigenvalue. 
However, the eigenvalue, jl2{T)^ is negative for all values of r G (—1, 1), and in order 
to obtain a bounded solution in the limit z — > 00 we set C2 = 0. The vector C, is the 
generalized left eigenvector satisfying A{t,Q)'^C, — (5(0)1 and is given by 

C= (-1/6,0,1/6)^. (3.77) 
At the boundary, the solution is 

ai(T,0) =col + ciC, (3.78) 
and the boundary condition (3.68) requires 

CO + y = 0, (3.79) 

so that ci — —bcQ. Thus, up to a single unknown constant, cq, which must be 
determined by matching, the boundary-layer solution is 

^bi(r,z) -co(l-6(C + zl)). (3.80) 

Because ^bi(T, z) is unbounded in the limit z — )■ 00, it is not possible to match it to the 
outer solution, ^out- We can think of the term, ^ -|- zl, in the boundary-layer solution 
is a truncated Taylor series expansion of the true solution around z = 0. To match 
the boundary-layer and outer solutions, a transition-layer solution is required for the 
strip of width 0(-\/e) along the boundary. 

Consider the stretched coordinate s = cr/^/e. Keeping terms to 0(e^/^), the 
transition-layer solution, ^tK^, s), satisfies 

VeGiO)^ + A{t, Of ^ti^O. (3.81) 
OS 

It is less clear how to truncate the above equation to obtain a leading order transition- 
layer solution. Because we must match the outer solution, 1, to the boundary layer 
solution that has the generalized eigenvector ^, we try a solution of the form 

^ti(T, s) = ao(T, s)l + ai(r, s)x{t, s), (3.82) 

where 

X{r,s)= ^ (l-T7i(T,Veg)), (3.83) 
v{t, y/es) 

and ao,i are unknown scalar functions. In the limit a — >■ 0, the deterministic flux 
across the boundary, 

>^{t, (j) = f{x{T, (T),y{T, cr)) - f{y{T, cr), x(t, a)), (3.84) 

(with f{x,y) given by (2.14)) vanishes and the eigenvector rji — 1. That is, 
the eigenvalue /Zi, corresponding to the eigenvector rji, vanishes on the boundary. 
Furthermore, it can be shown that 

limy(r,.) = -^^£||^C (3.85) 



Isolating intrinsic noise sources in a stochastic genetic switch 



17 



where 



da 



Substituting (3.82) into (3.81) yields 



y^G(O) 



ds 



ds 



da 

dx 

da 



aiAiT,Ofx = 0- 



(3.86) 



(3.87) 



To obtain the unknown functions ao.i(T, s) we project (3.87) with the right 
eigenvectors, ■0„(r, s), n = 0,1. After applying these projections (using the fact 
that iPqGx = 1, "^o-^ = 0, iI^qGI = v{t, a), and Gl = 0) and collecting leading 
order terms in e, we get 



dai ^ ^9ao 

— = v[T, Vesj^)- 
os OS 

— + s^i\ '{T)ai = 0, 
as 



(3.88) 
(3.89) 



where 



^'^ (r + 1) 



(3.90) 



-2 26 - 1) 

[t + W 

It turns out that p}-^^ is related to the curvature of the stability landscape normal to 
the separatrix; that is, if we define 

^{T,a) = ^{x{T,a),y{T,a)) (3.91) 

then 

/^^''^(^) = -£^(^'0)- (3-92) 

At = \/l ^ 26, the curvature vanishes and changes sign but is always negative 
(with fl'i{Tu) > 0) at the unstable fixed point, (t„,0). Divide the separatrix (ct = 
and — 1 < T < 1) into three regions: — 1 < r < 0, < t < Tq, and Tq < r < 1. 
The first region is ignored because it is in part of the domain V excluded from the 
stationary density function (see Sec. 3.1). The second region contains the unstable 
fixed point, and the third we can ignore as only extremely rare trajectories cross the 
separatrix in this region. Up to an unknown constant, the solutions to (3.88) and 
(3.89) are 



ao(T, s) 



exp 



1 



du, 



1 



ai(r, s) ^ a exp 
The transition layer solution is then 

/ M^^(r) 



lti('r,s) 



ei/2i>M(r) 7o 
1 



exp 



- exp 



x{t, s) 



dul 



(3.93) 
(3.94) 

(3.95) 
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The three solutions can now be matched. First, matching the transition layer 
solution to the boundary layer solution is done using the Van-Dyke rule. In terms of 
the boundary layer variable, z, the transition layer solution is 



Matching terms with the boundary layer solution yields 



a = bco 



t^Y\r) 



(3.96) 



(3.97) 



The composite boundary/transition layer solution is then 



lbl/tl(T, cr) = Co 



exp 



1 -■(<y)( \ 2 



dul (3.98) 



- exp 



-2^1 (r)- 



The final unknown constant, co, is determined by matching to the outer solution 
so that 

lim |bi/ti(r,s) = l, (3.99) 
which imphes that 

\/26Al"^(r) 

In order to evaluate the term in the numerator of the eigenvalue formula (3.13), we 
require the adjoint eigenfunction evaluated on the boundary (in a neighborhood of the 
unstable fixed point) which is 



Co = 



(3.100) 



6(T,0) 



2eA('^)(T) 



(-1 + 6C), re (G,\/l-26). (3.101) 



3.3. Principal eigenvalue 

We now have all of the components necessary to approximate the principal eigenvalue, 
using the formula (3.13). First, for the term in the denominator, we can approximate 
the adjoint eigenfunction with the outer solution ^o ^ 1 and the (unnormalized) 
stationary density with (3.15) (the higher order term ri can be ignored). Then 
the term in the denominator is simply the normalization factor, which can be 
approximated using Laplace's Method to get 



exp 



V 



^<^{x,y) - k{x,y) 
e 



27re 



ydrt(Z) 

where Z is the Hessian matrix (3.25) of $ at the stable fixed point {x^,,y 
that we have used the fact that fc(a;*,j/*) — ^(a;*,?/*) 



(3.102) 



Note 

and that the vector Vq is 



normalized so that its entries sum to one. 
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The term in the numerator of (3.15) requires the approximation (3.101) of 
the adjoint eigenfunction on the absorbing boundary. The integral can also be 
approximated using Laplace's Method, with 



(3.103) 



1 

71 



bV^~^2ef^[^\r) 



C^G(0)fo(r)exp 



-^$(t) - fc(T) 

e 



dr 



^-fe(x„,l/u) 



C'G(0)fo(T)exp 



e 



For convenience, we have defined functions on the boundary in the variable r with 

fo(r) = ro(x(r,0),2/(T,0)) (3.104) 

$(r) =$(x(T,0),y(r,0)) (3.105) 

fc(r) =fc(a;(T,0),y(r,0)), (3.106) 

where x{T,a) and y(T, cr) are defined in (3.62). 

Although the quantities ^{xu,yu) and k{xu,yu) must be computed numerically, 
the remaining unknown terms can be computed analytically by exploiting the reflection 
symmetry of the problem. Along F, we have that x = y and p — q so that the equation 
for $ and ro (3.16) can be written as 

i(T, 0) - A^r(r)F(T)] fo(T) = 0, (3.107) 
where we have defined 

^r(r) = -p(x(t, 0), 2/(r, 0)) = -q(x(t, 0), y(r, 0)). (3.108) 
The stability landscape function on the boundary is then 



l>(r) = $(x„,2/„) 



firiu)du. 



(3.109) 



The above is just an eigenvalue problem with three possible solutions, one of which 
can be excluded because there is a zero eigenvalue corresponding to the nuUspace of 
A. It can be shown [23] that if the diagonal elements of F{t) are such that at least two 
have oposite sign then only one of the remaining two eigenvalues has a corresponding 
positive eigenvector (fo(''') raust have positive elements) making the solution to (3.107) 
unique. It turns out that this is only true for r e (0, 1), which is due to the domain 
restriction caused by removing protein fluctuations. The result is 

+ 2t2 + 26(t - 1) + t 



Mr(r) = - - 

fo(r) = Q(l- 
We also have that 

|."(r) = -Mr) 



2t(1 



-T* + 2t^ + (3 + 26)r2 - 46r + 2b 



2t2(1-t)2 

And finally, using (3.77), (3.65), and (3.111) we get 



(3.110) 
(3.111) 

(3.112) 



C^G(0)ro(r) 



,(l-r). 



(3.113) 
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Combining these components together, we have the final result that 



1 



/2n 



(1 - T^)e-'=(="'"f") 



■ cxp 



(3.114) 



where t„ = a;„ + y„ — 1, and p,\ (r) is given by (3.90). As expected, the eigenvalue is 
exponentially small in e, which means that the height of the stability lanscape at the 
unstable fixed point, <i>(a;„, ?;„), must be approximated as accurately as possible. 

The remaining terms are often referred to as the 'prefactor', and except for the 
quantity k{xu,yu), slU of the terms in the prefactor (the derivatives of the stability 
landscape function and (3.113)) represent properties local to the fixed points. The 
remaining term in the prefactor depends on the function k{x,y), which depends 
on properties of the process not local to the fixed points and must be computed 
numerically. 



4. Results 



In this section, the results gathered throughout this paper are used to explore how 
removing the intrinsic noise that arises from protein production/degradation effects 
the random process. In particular, we examine the stability landscape and the 
metastable transition times. First, in Fig. 3 the numerical solutions to the ray 
equations (3.20) are used to generate level curves of the stability landscape function, 
$, for both the QD (3.18) and the full (3.60) Hamiltonians. For presentation, the level 




Figure 3. Level curves of the QD (blue) and full process (orange) stability 
landscape function <I>. Black lines show the domain boundary for the QD process, 
and the dashed line shows the boundary of positive protein levels, (a) The QD 
result is compared to the full process with </? = 0.1 so that the protein noise is 
weak compared to promoter noise, (b) Same as (a) but with ip = 1 so that the 
protein noise strength is comparable to promoter noise. Parameter values are the 
same as Fig. 2. 

curves are shown in the {x + y,x — y) plane, with the separatrix along the left edge 
{x — y = 0) of each frame. In Fig. 3a, level curves for the full process are shown for 
= 0.1 so that the protein noise is small compared to promotor noise. Recall that the 
parameter ip controls the strength of protein noise relative to the strength of promotor 
noise so that there is no protein noise in the limit (/3 — ^ 0. The resulting curves match 
closely in a neighborhood of the stable fixed point and extend out toward the unstable 
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saddle point. As expected, the level curves begin to diverge the farther away from 
either fixed point they are. In Fig. 3b, the strength of the protein noise is increased, 
with if = I, and in this case, the stability landscape of the two processes are quite 
different, matching only near the fixed points of the deterministic dynamics. 

As a result of the domain restriction effect, level curves of the full process extend 
into regions the QD process is excluded from. This is because protein levels can only 
cross above the line y — I — x, not below it, when there is no protein noise. This is 
also true of the lines x = 1 and y = 1. The domain restriction effect is eliminated 
when protein noise is added back into the process, even if it is very small compared to 
promotor noise. For the model considered here, the effect is of no serious consequence 
for metastable transitions as the restricted domain still contains all three fixed points. 
However, a model of a more complex gene circuit might be significantly affected by 
removing protein noise — especially if this restricts the domain for the protein levels 
in such a way as to generate qualitatively different behavior, which would imply a 
nontrivial contribution of protein noise, no matter how negligible it may be. 

Because of the symmetry in the problem, we obtained analytical results for various 
quantities on the separatrix, including the shape stability landscape. We can use these 
results to obtain an analytical approximation of the probability density for the position 
along the separatrix a trajectory passes through as it transitions from one basin of 
attraction to another. Using the results of Sec. 3.3, the stationary density along the 
separatrix is given by 



fo(T)e exp 


-i$(r)' 






i|.(r) 


dr 



Pexit(T) - — — (4.1) 

e-'=(^) exp -^$(r) dr 

where we remind the reader that t = x + y — 1 and x = y along the separatrix. The 
only term that cannot be obtained analytically is the function A:(t), which can be 
ignored as a first approximation. For simplicity, we also average over the promotor 
state to get the scalar marginal probability density for the exit point. Then, using 
Laplace's method, the exit density is 



^oxit(T)^ W^!^exp 



27re 



i(|.(r)-$(r„)) 



re (0,1), (4.2) 



where <1>(t) is given by 



|.(r) = - ^ - 61og(T) - 2 log(l ~ r). (4.3) 

The QD exit density approximation is shown in Fig. 4 along with two histograms 
obtained from Monte-Carlo simulations. While the histogram for if = 0.1 is close to 
the QD approximation, it is evident that some trajectories pass through the separatrix 
in the interval (—1,0), which is impossible without protein noise. This effect becomes 
negligible when the protein noise is reduced to ip = 0.01. 

The FETD for a trajectory, starting from a stable fixed point, to reach the 
separatrix is asymptotically exponential in the large time limit, and the timescale 
is determined by the principal eigenvalue, Aq, from equation (3.114). We can then 
approximate the mean exit time with T ~ I/Aq. In Fig. 5 the mean exit time is shown 
on a log scale as a function of 1/e along with results from Monte-Carlo simulations. 
Notice that the approximation and the Monte-Carlo simulations are asymptotically 
linear as 1/e — >■ oo. For the approximation, the slope of this line is determined by the 
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Figure 4. The density of exit points along the separatrix as a function of 
r = x + y — l, with x = y. The black curve shows the analytical approximation for 
the QD {(p = 0) process, and the symbols show histograms from 10* Monte-Carlo 
simulations for different values of (p. 




Figure 5. The mean exit time on a log scale as a function of Symbols 
represent 10^ averaged Monte Carlo simulations for various values of (p. The 
dashed black line is the approximation, T = 1/Ao, for the QD process (i^ = 0). 



height of the stabihty barrier, while the prefactor affects the vertical shift. 

From the Monte Carlo results (symbols with grey lines) we see that the mean exit 
time converges to the approximation as — >■ 0. However, it is clear that the slope 
of the analytical curve is slightly different then that of the Monte Carlo results even 
when (fi is small. Thus, we may think of the mean exit time approximation for the 
QD process as an asymptotic approximation of the full process in terms of the small 
parameter <C 1 so long as e is also small but not too small. 

5. Discussion 

Understanding how different noise sources affect the dynamics of a gene circuit is 
essential to understand how different regulatory components interact to produce the 
complex variety of environmental responses and behaviors. Even if one excludes 
extrinsic noise sources — such as environmental and organism-to-organism variations — 



Isolating intrinsic noise sources in a stochastic genetic switch 



23 



there are several sources of intrinsic noise, such as fluctuations in translation, 
transcription, and the conformational state of DNA regulatory units. 

The behavior we are interested in understanding is a transition from one 
metastable state to another. Each metastable state corresponds to the stable steady- 
state solutions of the underlying deterministic system. The bistable mutual repressor 
model has two identical stable steady states separated by an unstable saddle node. On 
small time scales, the protein levels fluctuate near one of the two stable steady states. 
On large time scales, fluctuations cause a metastable transition to occur, where the 
protein levels shift to the other steady state by crossing the separatrix containing the 
unstable saddle point. 

To understand how metastable transitions can be induced by promotor noise, we 
consider a discrete stochastic model of a mutual repressor circuit where the protein 
levels change deterministically, which we call the QD process. This is done by fixing 
the promotor state and then taking the thermodynamic, or large system size, limit of 
the protein production/degradation reactions. We then compare the QD process to 
the full process that includes protein noise. 

We find important qualitative differences that persist even when the magnitude of 
the protein noise is small compared to promotor noise. In particular, without protein 
noise and after initial transients, the protein levels are restricted such that the total 
amount of protein is never less than half its maximum possible value; that is, the total 
number of proteins is such that n + m > a/S, where a is the protein production rate 
and S is the degradation rate. Said another way, assuming that the random process 
starts at one of the deterministic stable fixed points, promotor fluctuations could never 
push the protein copy numbers so that, for example, only a single copy of each protein 
is present. In contrast, the stability landscape for the full process and Monte-Carlo 
simulations show that protein levels are able to reach all positive values. While this 
restriction does not stop the QD process from exhibiting metastable transitions, more 
complex gene circuits may require protein fluctuations, even if they are very small, in 
order to function correctly. 

Appendix 

In this appendix, we summarize the numerical methods and tools used throughout 
the paper. Most numerical work is performed in Python, using the Numpy/Scipy 
package. For more computationally-expensive tasks, we use Scipy's Weave package to 
include functions written in C, which allows us to use the GNU Scientific Library for 
numerical Integration of the ray equations (3.20), and for random number generators 
used in Monte Carlo simulations. 

There are a few notable observations regarding integration of the ray equations. 
First, characteristic projections, {x{t),y{t)), have a tendency to 'stick' together along 
certain trajectories, peeling off one at a time (see Fig. 2). To adequately cover the 
domain with rays, a shooting method must be used to select points on the Cauchy 
data. For more details on this see Ref. [13]. We found that the simplest method was to 
use the secant method (we use the "brentq" function in the Scipy.integrate package) 
to minimize the euclidian distance between the final value of {x{t),y{t)) along the 
separatrix and the saddle node. This method is convenient since it does not require 
knowledge of the Hessian matrix, Z{x,y). Second, the value of w used to generate 
Cauchy data must be chosen small enough to get accurate results. However, we found 
that if it is chosen too small, rays are no longer able to cover the domain, and more 
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importantly, we could no longer generate a ray that reaches the unstable fixed point 
on the separatrix. 

For the mutual repressor model, the trajectory connecting the stable fixed point 
to the saddle is one of the curves along which characteristics tend to stick to each 
other. Suppose that 9m is the point on the Cauchy data (3.24) that generates the 
ray that connects the fixed points. Then small perturbations = 9m + 59, \59\ <C 1, 
cause the characteristic, {x{t), 2/(0)' to diverge sharply away from the saddle. This not 
only makes it difficult to compute 9m, but also creates difficulties for computing the 
function k{x,y) (see Sec. 3.1.1). Since the expanded set of ray equations (3.50) track 
the derivatives of x, y, p, and q with respect to the point on the Cauchy data, which, 
for values of 9 near 9m, becomes very large as the ray approaches the saddle point. 
As the expanded variables become very large, computing Z using equation (3.49) is 
unstable. Furthermore, this effect becomes worse as the initial value, uj = $(0), goes 
to zero. 
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